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Abstract — A nonperturbative analytic solution of the high charge and energy (HZE) Green s function is 
used to implement a computer code for laboratory ion beam transport in multiple-layered materials. The 
code is established to operate on the Langley nuclear fragmentation model used in space engineering 
applications. Computational procedures are established to generate linear energy transfer (LET) distri- 
butions for a specified ion beam and target for comparison with experimental measurement. Comparison 
with 56 Fe ion with Pb-Al and Pb-(CH 2 ) T targets shows reasonable agreement. 


1. INTRODUCTION 

Green’s functions were identified as the likely means 
of generating efficient HZE shielding codes for space 
engineering which are capable of being validated in 
laboratory experiments (Wilson ei a/., 1989). A deri- 
vation of the Green’s function as a perturbation series 
gave promise for development of a laboratory-vali- 
dated engineering code (Wilson et al. , 1990) but 
computational inefficiency provided a major obstacle 
to code development (Wilson and Badavi, 1992). 
More recently, nonperturbative approximations to 
HZE Green’s functions have shown promise in 
providing an efficient validated engineering code 
(Wilson et al. , 1993a, c). Previous work has found a 
solution to HZE transport in a homogeneous 
medium using nonperturbative methods (Wilson and 
Badavi, 1992; Wilson et al. , 1993b, c). In the present 
report, we derive solutions for inhomogeneous multi- 
layered media. The resulting computer code is used to 
derive LET spectra behind multilayered targets for 
ion beams with Z ^ 28 corresponding to the major 
components of the galactic cosmic ray spectrum. The 
results of the computation are compared with 56 Fe 
accelerator beam experiments with Pb-Al and 
Pbi-(CH 2 ) T shield configurations. 

2. GREEN’S FUNCTION FOR A SINGLE 
MEDIUM 

We restrict our attention to the multiple charged 
ions for which the Boltzmann equation may be 
reduced (Wilson, 1977a) to: 

f - 4 $(£) + <x,1 <t>,(x, E)-I a jk <P k (x, £), 
ex cE J * 


where ^(x, £) is the ion flux at x with energy £ 
(MeV, amu), $,-(£) is the change in E per unit dis- 
tance, a } the total macroscopic reaction cross section 
and a* the macroscopic cross section for collision of 
ion tvpe k to produce an ion of type j. The solution 
to equation (1) is to be found subject to the boundary 
condition: 

*;(0 ,£)=#£). (2) 

w'hich for laboratory beams has only one value of j 
for which fj(E) is not zero and that /•(£) is described 
by a mean energy £ 0 and energy spread a such that: 

f.(E) = —^L- exp[ — (£ — £o) J /2ffT O) 

V 2na 

The usual method of solution is to proceed solving 
equation (1) as a perturbation series (Wilson 
19 a. b; Wilson et al., 1990). In practice, the compu- 
tational requirements limit the usefulness of the tech- 
nique for deep penetration (Wilson and Badavi, 
19921 . 

Tee Green’s function is introduced as a solution of: 

[~2_ _ ^ S,{E) + <r y J G jm {x, E, £„) 

= I <r jk G km (x, £.£,). (4) 

k 

suNcd to the boundary condition 

G ym (0,£,£ o ) = M (£-£„)• ( 5 > 

Tee solution to equation (1) is given by superposition 
as 

*,(*,£)- 1 Jc yt (.x, £,£')/*(£') d£'. (6) 


mi 


0 ) 
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If G lK (x, £, £') is known as a transcendental func- 
tion, the evaluation of equation (6) may be accom- 
plished by simple integration techniques, and the 
associated errors in numerically solving equation (1) 
are avoided (Wilson et al., 1991). 

The above equations can be simplified by trans- 
forming the energy into the residual range as: 

rj-j' dE'lSjE'), 

and defining new field functions as: 

fj) = Sj(E)<t>j(x, E) 

#,„(*, r Jt r' m ) = ■?,(£)(?, m (x, £, £') 

fM) = Sj(E)fj(E) 

and equation (4) becomes: 


In the above, the g -function of n -arguments is found 
recursively by: 

g(j) = (19) 

and 

g(.j\*Jl' Jir>jn + l) 

g (J { , Jl jn — 1 i jn ) g^.j\ » Jt jn— ].jn+\'} 


(7) 

( 8 ) 
(9) 

( 10 ) 




( 20 ) 






dx dr 




= 1 -^k„(x,r k ,r’ m ), (11) 

k V k 

where i) is the range scale factor as v/j = v m r m and is 
taken as v. = Z~ jA t and the boundary condition is 
now: 

3>(0, r r r' m ) = 5 jm 6( rj - r' m ) (12) We note that: 

and with solution to the ion fields given by 

'l> j (.x,r J ) = Y 


Note that the #£(*,/>, r' m ) are purely dependent 
on x for /> 0 which we represent as ^(*) 
(Wilson and Badavi, 1992). In terms of the above, the 
solution to equation (1) becomes (Wilson and Badavi, 
1992) 

^x,r J ) = t-’^f J (r J + x) 

+ M(*)[^(C)-A(r;j]. (21) 

In equation (21), and r’ m , are given by the upper 
and lower limits of the inequality (17). The symbol 
f m (r' m ) refers to the integral spectrum: 


fjiriiX’ r j* dr’ m . (13) 


with 


The solution to equation (11) is written as a pertur- 
bation series: 


P m ( r m) — j */*.( r ) dr. 
Pjr' m ) = F m {E’) 

-J> 


F m (E') 


(£) d£ 


and 


Vjmix, r r r „) = £ #£(*, G, r'J 


(14) 


w'here 


h 


d£/S m (£). 


( 22 ) 


(23) 


(24) 


(25) 


&j°J(x, fj, r' m ) = g(j ) 6 jm 6(x + Tj — r' m ) (15) 


and 


(16) 


(17) 


where ^"(x, r } , r' m ) is zero unless 

P, , V. 

— (r ; + X) <£ r m «£ rj + X 
t’m 

for v m >Vj. If i)> v m , as can happen in neutron 
removal, the negative of equation (16) is used and the 
upper and lower limits of equation (17) are switched. 
The higher terms are approximated as: 

&»(x,r r r'J* £ 

Jkj >3 — - 1 


We now introduce nonperturbative terms for the 
summation in equation (21). 

First, we recall that the g-function of n -arguments 
was generated by the perturbation solution of the 
transport equation neglecting ionization energy loss 
(Wilson et al 1989) given by: 

(26) 

subject to the boundary condition: 

gj*(0) = 6 Jn , (27) 

for which the solution is 

g J „(x) = S jm g(m) + a jm g(j, m) + -- (28) 

It is also true that: 

g,n, (* ) = I ft* (* - y >ft" O’ ) (29) 


x(v m - V,) 


(18) for any positive values of x and y. Equation (29) may 
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be used to propagate the function g im (x) over the 
solution space, after which: 


$,Jx, r r rj ft e - *'* S, m 5(x + r, - rj 

+ v,\g )m (x) - e-*' J T 5 jm }/x{v„ - v,). (30) 


The approximate solution of equation (1) is then 
given by 


+ v,\g, m (x) - e-y _M 


(31) 


which is a relatively simple quantity (Wilson et al., 
1993a). 


3. GREEN’S FUNCTION IN A SHIELDED 
MEDIUM 

The major simplification in the Green’s function 
method results from the fact that the scaled spectral 
distribution of secondary ions to a first approxi- 
mation depends only on the depth of penetration as 
seen in equations (16), (18) and (30). Our first ap- 
proach to a multilayered Green’s function will rely on 
this observation and assume its validity for multilay- 
ered shields. 

Consider a domain labeled as 1 which is shielded 
by a second domain labeled as 2; the number of type 
j ions at depth .v in 1 due to type m ions incident on 
domain 2 of thickness y is: 

(32) 

* 


*«*(*.* T P r m) ~ e-''' T -" ! ' v 5 jm 

X 5(x + pv + r ] - rj + v) 

_ c -»i/ T -« r 2 /> 3 im ]j(x + py)(v m - v t ). (37) 

Equation (37) is our first approximation to the 
Green’s function in a shielded medium (two layers) 
and is easily modified to multiple layers (see Appen- 
dix). We now consider the first spectral modification. 

It is easy to show that the first collision term has 
the properties: 


y\'lJx, y, r,, rj 


v>o 


7 iy<" 


»-®l 




for r' = r ' 


= f , , 
J - :t for r _ = r 


I”* ~ v j\ 


(38) 


We use these properties to correct the average spec- 
trum as: 


&ilJx, V, r,, r’J = 


V,gl2,m(X,y) 

\v m - r,|(.v 4- pv) 

+ b jm (x,y)(r’ m -fJ , (39) 


where g\ n Jm (x, y) is the first collision term of equation 
(37) and* 


r m (/ mu 4 " r ml)l- 


(40) 


is the midpoint of f’ m between its limits given by 
equation (36). The b Jm term of equation (39) has the 
property that: 


i: 


bj„(x,y)(r' -rj dr’ = 0 , 


(41) 


The leading term in equation (32) is the penetrating 
primaries as: 

gv. im {x,y) = ^-^b jm 

+ \gr. im (x,y)-t-^-^5 im ), (33) 

where all higher order terms are in the bracket of 
equation (33). 

The first term of the scaled Green’s function is 
then: 

X 5[x + T } - {r' m - py)], (34) 

where p is the range factor for the two media: 

p = Ri, (E)IR 2J (E). (35) 

The ratio is shown for protons in Fig. 1. We take a 
single value for p corresponding to 600MeV/amu. 
The secondary’ contribution is similarly found by 
noting that equation (17) becomes: 

— (r 4- x 4* py ) ^ r' m < — r } + x + py , (36) 

v m v * 


ensuring that the first term of equation (39) is indeed 
the average spectrum as required. The spectral slope 
parameter is found to be: 

b Jm (x,y) = y ; u m (<r lym e- tf, "' T " ff2 «- r - a ljm 

xt~*u*—iiy)i[(x + py){v m - v ; )\v m - (42) 

A similarly simple spectral correction could be made 
to the higher order terms. The spectral correction 
given in equation (42) is included in the present 
Green’s function code. 



from which the average spectrum is evaluated. The 
full approximate Green’s function is then: 


Fig. 1. Ratio of range in water to range in aluminum for 
proton beams. 
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4. LET SPECTRA FOR LABORATORY BEAMS 

We use the boundary condition appropriate for 
laboratory beams given by equation (3). The cumu- 
lative spectrum is given by: 




1 -erf 


E-E o 


(43) 


/J 

The cumulative energy moment needed to evaluate 
the spectra! correction is: 


Zj(E)= ] -E o 


1 — erf 


E-E a 




a /J 


+ 




exp 


(£-£ 0 ) : 
2c 2 


. (44) 


The average energy on any subinterval (£,,£;) is 
then: 


E - [£,(£,)- £,(£; )]/[£/£])- £>(£:)]- (45) 

The beam generated flux is: 

«/0(*> >'< ',) = e-«'-*"/(r, + * + Pj) 

m./ 

+ X (* X » ~ ^*1 

m 

x[/L(rl)-/ ( 46 ) 

where £ is evaluated using equation (45) with £, , and 
as the lower and upper limits associated with r m} 
and r' mu . 

A series of evaluations for a lead scattering foil 
(2.24 g/cm 2 ) in front of a water target is shown in 
Fig. 2. The lead scattering foil is usually part of the 
accelerator beam line so that the fragments from the 
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Fig. 2 (a) and (b). 
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Fig. 2. Differential fluence for 525 MeV/amu *Fe beam with a 2.5 MeV/amu standard deviation after 
passing through a 2.24 g; cm 2 lead scattering foil and a water target, (a) Ocm H 2 0; (b) 5 cm H : 0; 

(c) 10 cm H 2 0. 


lead target are seen as contamination. Clearly, these 
fragments must be modeled to properly interpret the 
attenuation of the beam in the water target in actual 
experiments. 

5. NUCLEAR DATA BASE 

The nuclear absorption cross sections are fits to 
quantum mechanical calculations developed at the 
Langley Research Center over the past 20 years 
(Wilson, 1973, 1974; Wilson and Costner, 1975; 
Wilson and Townsend, 1981; Townsend and Wilson, 
1986) and are considered reliable to about 10%. The 
nuclear fragmentation cross sections for most nuclei 
on hydrogen targets are taken from Silberberg et al. 
(1983) and are augmented for light fragment pro- 
duction with the Bertini model (Bertini, 1969). It 
was noted that early versions of these cross sections 
failed to conserve mass and charge (Wilson et a/., 


x 10" 1 
l.o r- 


^ 0.8 \- 
§ 0.6 


b 

« 0.4 
C 

m 

Cl. 


0.2 


O (IA j0jp )/A p 


□ a 




Ejn°° 


" n0 


. 6 °° 




_L 


_L 


4 8 12 16 20 24 28 

Projectile ion charge 


Fig. 3. The absorption cross section in hydrogen targets 
and mass averaged production cross sections at 
600 MeV/amu for various projectiles. 


1974) and still exhibit mass loss for 10 ^ Z ^ 22 by 
as much as 30%. This is displayed in Fig. 3 where 
< 7 abs is compared to TLA,a (p !A p , where A,- is fragment 
mass, c,p is the fragmentation cross section for 
projectile p and A p is the projectile mass. The 
breakup of light nuclei ( A < 4) is taken from the 
quantum calculations of Cucinotta et al. (1993). The 
fragmentation of the remaining nuclei ( A p > 4) is 
evaluated from the latest versions of the 
NUCFRAG model (Wilson et al. 1987a, b). Since 
the public release version of NUCFRAG (HZE- 
FRG1, Townsend et al ., 1993), a de-excitation 
scheme for mass two and mass three fragments and 
a coulomb trajectory calculation have been added 
for more realistic cross sections at low energy 
(Wilson et al 1993a), The elemental fragmentation 
cross sections are displayed in Fig. 4 at several 
energies. The reduced light fragment production at 
low energy results from coulomb trajectory correc- 
tions. This is the same data base used in the most 
recent energy dependent engineering code HZETRN 
(Shinn et al 1992). 

The transport codes usually represent a reduced 
set of isotopes. In the past, we usually represented 
each charge group with an associated mass taken as 
the nearest mass on the stability curve for the given 
fragment charge. The most recent version of 
HZETRN uses an isobaric flux representation with 
the nearest charge on the stability curve and the 
distance to the nearest isobar was calculated 

D = (A t - A,) 2 + 4(Z, - Z { )\ (47) 

where A iy Z, represent the fragment and A h Z, 
represent the listed isobar mass used in the calcu- 
lation and nearest charge to the stability curve. The 
present calculation uses an 80-isotope representation 
and the nearest isotope in the list is found using 
equation (47). 
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FiG. 4. The fragment production cross section in H 2 0 targets including coulomb corrections. 


6. EXPERIMENTAL METHODS AND 
COMPARISON 

The 56 Fe nuclei were accelerated to 600 MeV/amu 
at the Lawrence Berkeley Laboratory Bevalac facility 
and passed through a series of beam transport 



Fig. 5. Measured CR-39 response parameter (Benton et al , 
1986). 


elements, triggering devices and a 2.24 g/cnr lead foil 
prior to exiting the beam pipe and impacting the 
target. Two targets of 2 g/cm 2 Al and 4.6 g/cm 2 of 
polyethylene (CHj)* were used to evaluate their 
transport properties. The beam energy is inferred to 
be 557 MeV/amu when only the lead foil and target 
are considered for transport analysis. The trans- 
ported beam exciting the target was measured using 
CR-39 plastic foils (Benton ei al . , 1986). The beam 
intensity was measured by a monitoring foil in front 
of the target. The detectors and targets are run in 
good geometry so that acceptance corrections are not 
required. 

The detector response is assumed to be approxi- 
mately Gaussian with an LET dependent width T 
shown in Fig. 5. A correction for non-Gaussian 
contributions is taken as: 

R(L, Lo) = 0.8 J— t ~' L - 

s/271 <75 

+ 0.2-d=e- ,i - w * 3, t (48) 

s/2 no] 

where <r 0 = 0.4247 f and <r, (taken as 2.4ff 0 ) is fit to the 
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Fig. 6. Transport experiments for a lead-shield multilayer combination for 557 MeV/amu iron beams, 
(a) 2.24 g/cm 2 Pb + 2 g/cm 2 Al; (b) 2.24 g/cm 2 Pb -b 4.6 g/cm (CH 2 ) T . 


high LET side of the primary ion peak. The response 
function of equation (48) is used to compare the 
theory to the experiment. 

The distribution of ions produced in passing a 
557 MeV/amu *Fe beam through a 2.24 g/cm 2 of lead 
and the two target materials (separately) was mapped 
into detector response using equation (48). The com- 
parison with experimental measurements is shown in 
Fig. 6. While the calculated result for polyethylene is 


in good agreement with the experimental data 
(Fig. 6(b)), the calculated aluminum curve (Fig. 6(a)) 
suggests that the aluminum fragmentation cross sec- 
tions may be 20-30% low. 
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APPENDIX 

The preceding formalism is extended to a three-layer 
configuration as follows. The solution to equation (26) in a 
three-layered medium is: 

gm^(x,y, :) = 'Zg i j t (.x)g nl tY)g ilm (:). (A!) 
*/ 

The leading term is the penetrating primaries, and equation 
(Al) may be written as: 

+ y. -) - <5jJ. (A2) 

The scaled Green’s function is then: 

9W*. * V r -> * — S„ 

*<5(* + Pj.v + p 3 r + r- r'J + v y [ ?l j jjm (x. y, :) 

- e —'*-*"-* 1 S^Kx + Pl y + p 3 --)(r„ - 1,), 

(A3) 

where p. = R xj {E)jR v (E) and p 3 = A,j(£)/A v (£). The 
range condition of equation (17) becomes: 

^-(r + x + p 2 y + PjZ) < r' m s£ — r ■ + x + p : + p,r. (A4) 

I'm 

The spectral corrections are similarly derived. 
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